In 2012, we placed 96 dataloggers recording temperature and humidity across Table Mt. and Silvermine. The loggers were stratified across gradients of elevation, slope, exposure and hillslope position (hilltop, valleybottom), and distance from both False Bay and the Atlantic. In previous analyses, the data has been analyzed to extract diurnal Tmin, Tmax, VPmax (maximum daily vapor pressure), and RHsat.hrs (hours per day with relative humidity above 95%).
All of these analyses and the summary data are available at: https://github.com/dackerly/table-mt-microclimate-2012.git
General summary statistics about the sites and daily weather values can be found in ‘manuscript-summary-stats.r’.
The goal of these analyses is to examine spatial variation in daily weather values, running models for each day through 2012. The initial question is whether the important spatial variables vary with season, but as shown below a key result is enormous day to day variation in the ability to explain spatial variation in weather variables. I’m still searching for explanations in regional weather or other factors to understand why the models differ so much from day to day.
#----------------#
# LOAD LIBRARIES #
#----------------#
rm(list=ls())
library(raster)
## Loading required package: sp
library(RColorBrewer)
#----------------#
# LOAD DATA #
#----------------#
‘meta’ contains metadata on the 96 sites. ‘topo10’ and ‘topo30’ have a series of topographic features sampled from 10m and 30m dem rasters, respectively, based on lat/long positions of each site. SiteID is the unique name for each site.
Some of the variables in ‘meta’: * elevation, slope, and aspect: measured in the field - these are generally not used, and instead we use a consistent set of values taken from the 10m dem. * use4Analyses, use4ClimateSummaries: sites excluded due to incomplete data (see ms for criteria)
Variables in topo10 and topo30: * d2at, d2fb, d2cs: distance to atlantic, distance to false bay, distance to coast = min(d2at,d2fb) * rad080, rad172, rad355: solar radiation at summer and winter solstice, and equinox. Not generally used as we have diurnal radiation at each site * thl0, thl315, thl337: Topographic heat load as defined by McCune and Keon 2002 JVS, with the axis for calculation of the ‘folded aspect’ (and thus the maximum heat load) at 0, 337 or 315 degrees. The logic is that in the S Hem, NW facing slopes would be hotter than N facing, due to higher afternoon temperatures. * plow050, 125, 250, 500: Percent lower pixels - the percent of pixels surrounding a focal pixel with a lower elevation, in radii of 50, 125, 250 or 500 m. A measure of potential for cold air to flow away from a spot. Higher values indicate hilltops * tpi050, 125, 250, 500: Elevation of a focal pixel - mean elevation in the specified radius. Another measure of hillslope position, with positive values for hilltops, negative for valley bottoms.
meta <- read.csv('data/csv_masters/location_meta.csv',as.is=T)
dim(meta)
## [1] 96 27
head(meta)
## domain siteID logger UTM.east UTM.north LON_DD LAT_DD
## 1 SILVERMINE CHK1 10008095 256907 6228779 18.36633 -34.05353
## 2 SILVERMINE CHK2 10006349 257070 6229153 18.36820 -34.05020
## 3 SILVERMINE CHK3 10008094 257271 6228903 18.37031 -34.05250
## 4 SILVERMINE CHK4 10008096 257420 6229224 18.37201 -34.04964
## 5 SILVERMINE CHK5 10008093 258353 6229582 18.38221 -34.04663
## 6 SILVERMINE CHK6 10006348 258397 6229772 18.38274 -34.04493
## elevation slope_deg aspect_deg magn.true.aspect use4Analyses
## 1 120 19 265 M 1
## 2 326 8 350 M 1
## 3 325 32 240 M 1
## 4 367 26 355 M 1
## 5 436 35 15 M 1
## 6 455 33 280 M 1
## use4ClimateSummaries route CollSeq series deployed.as deployXssn
## 1 1 CB 1 83 10008095 0
## 2 1 CB 2 84 10006349 0
## 3 1 CB 3 85 10008094 0
## 4 1 CB 4 86 10008096 0
## 5 1 CB 5 87 10008093 0
## 6 1 CB 6 88 10006348 0
## start.date start.hour end.date end.hour log.max.d VegSurveyPlot
## 1 11/10/11 18 12/31/13 23 150 111110-A-09
## 2 11/10/11 18 12/31/13 23 150 111110-A-07
## 3 11/10/11 18 12/31/13 23 150 111110-A-08
## 4 11/10/11 18 12/31/13 23 150 111110-A-06
## 5 11/10/11 18 12/31/13 23 150 111110-A-04
## 6 11/10/11 18 12/31/13 23 150 111110-A-03
## VegSurveyDate distOnTransect Orig.siteID
## 1 11/10/11 NA
## 2 11/10/11 NA
## 3 11/10/11 NA
## 4 11/10/11 NA
## 5 11/10/11 NA
## 6 11/10/11 NA
topo10 <- read.csv('data/csv_masters/topo10.csv',as.is=T)
dim(topo10)
## [1] 96 22
names(topo10)
## [1] "siteID" "domain" "elevation" "slope" "aspect"
## [6] "d2at" "d2fb" "d2cs" "rad080" "rad172"
## [11] "rad355" "thl0" "thl315" "thl337" "plow050"
## [16] "plow125" "plow250" "plow500" "tpi050" "tpi125"
## [21] "tpi250" "tpi500"
topo30 <- read.csv('data/csv_masters/topo30.csv',as.is=T)
dim(topo30)
## [1] 96 21
names(topo30)
## [1] "siteID" "domain" "elevation" "slope" "d2at"
## [6] "d2fb" "d2cs" "rad080" "rad172" "rad355"
## [11] "thl0" "thl315" "thl337" "plow050" "plow125"
## [16] "plow250" "plow500" "tpi050" "tpi125" "tpi250"
## [21] "tpi500"
‘dw’ has all daily weather values for all sites and days. The dlySummary contains summaries of weather data on each day - average weather across sites. dlySol has diurnal radiation at each site for each day of the year, from Adam Wilson’s analyses on a 30m dem.
dw <- read.csv('data/csv_masters/2012daily.csv',as.is=T)
head(dw)
## doy11 year month day Tmin Tmin.time Tmax Tmax.time Tmean RHmin
## 1 366 2012 1 1 14.936 6.166667 25.210 17.16667 20.0730 71.361
## 2 366 2012 1 1 13.762 6.166667 27.161 15.00000 20.4615 78.299
## 3 366 2012 1 1 13.233 6.000000 22.776 17.33333 18.0045 83.567
## 4 366 2012 1 1 13.546 6.000000 25.647 15.33333 19.5965 81.918
## 5 366 2012 1 1 12.413 6.000000 22.896 13.83333 17.6545 91.167
## 6 366 2012 1 1 12.727 6.000000 27.014 14.00000 19.8705 84.421
## RHmax VPmax RHsat.hrs date siteID doy
## 1 77.590 2.239716 0 2012-1-1 CHK1 1
## 2 82.556 2.267693 0 2012-1-1 CHK2 1
## 3 85.799 2.215489 0 2012-1-1 CHK3 1
## 4 83.437 2.133584 0 2012-1-1 CHK4 1
## 5 87.839 2.011852 0 2012-1-1 CHK5 1
## 6 87.179 2.199438 0 2012-1-1 CHK6 1
dlySummary <- read.csv('data/csv_outfiles/dlySummary.csv',as.is=T)
head(dlySummary)
## doy year month day date Tmin Tmax Tmean dtRange
## 1 1 2012 1 1 2012-01-01 11.39483 21.62492 16.50988 10.230091
## 2 2 2012 1 2 2012-01-02 14.97574 19.01598 16.99586 4.040239
## 3 3 2012 1 3 2012-01-03 14.75741 23.10118 18.92930 8.343773
## 4 4 2012 1 4 2012-01-04 10.78530 23.34831 17.06680 12.563011
## 5 5 2012 1 5 2012-01-05 13.78350 25.12926 19.45638 11.345761
## 6 6 2012 1 6 2012-01-06 14.36026 27.59257 20.97641 13.232307
## RHmin RHmax VPmax RHsat.hrs sommax sommin Kps.mean_hPa Kppt_mm
## 1 83.64332 92.68127 1.920619 5.545455 16 15 997.0 0.6
## 2 52.31199 99.71526 1.937674 19.000000 15 14 995.5 11.6
## 3 52.63628 98.78947 1.702166 15.679924 7 7 997.2 0.2
## 4 39.66210 94.06428 1.792897 4.657197 11 11 998.2 0.0
## 5 28.20923 95.26178 1.872311 7.526515 11 7 999.6 0.0
## 6 51.95174 85.88167 2.043691 1.242424 11 8 999.4 0.0
## Kwspd.mean_m.s Kwspd.max_m.s
## 1 2.991667 5.6
## 2 2.504167 5.4
## 3 2.279167 4.7
## 4 2.304167 3.1
## 5 1.975000 2.7
## 6 2.170833 4.8
dlySol <- readRDS('data/Rdata/dlySol.Rdata')
dlySol[1:6,1:6]
## rad_tot_001 rad_tot_002 rad_tot_003 rad_tot_004 rad_tot_005
## CHK1 8663.728 8652.467 8640.090 8626.597 8611.990
## CHK2 9244.997 9240.468 9235.409 9229.817 9223.687
## CHK3 7737.912 7724.081 7708.918 7692.427 7674.611
## CHK4 8369.773 8371.854 8374.001 8376.205 8378.459
## CHK5 8767.362 8767.380 8767.277 8767.045 8766.678
## CHK6 8059.166 8050.924 8041.838 8031.906 8021.127
## rad_tot_006
## CHK1 8596.270
## CHK2 9217.014
## CHK3 7629.731
## CHK4 8380.751
## CHK5 8766.166
## CHK6 8009.500
Before fitting any spatial models of weather variation, we want to check for correlations and spatial autocorrelation of predictors. Based on the topo10 variables, here are all variables with pairwise correlations R2 > 0.5. Commented code can be run to see more detail. As expected, radiation load variables are mostly correlated, and hillslope position variables are correlated. Below, we run alternative models swapping these in and out, but not using them in the same model.
# cor(topo10[,-c(1:2)],use='pair')
# pairs(topo10[,-c(1:2)])
# names(topo10)
#
## find correlations > 0.5
i <- 3; j <- 4;
for (i in 3:(ncol(topo10)-1)) for (j in (i+1):ncol(topo10)) {
rr <- cor(topo10[,i],topo10[,j],use='pair')
if (abs(rr)>sqrt(0.5)) print(c(names(topo10)[c(i,j)],round(rr,3)))
}
## [1] "slope" "rad355" "-0.746"
## [1] "rad080" "rad172" "0.952"
## [1] "rad080" "thl0" "0.827"
## [1] "rad080" "thl337" "0.78"
## [1] "rad172" "thl0" "0.723"
## [1] "thl0" "thl315" "0.803"
## [1] "thl0" "thl337" "0.951"
## [1] "thl315" "thl337" "0.946"
## [1] "plow050" "plow125" "0.814"
## [1] "plow125" "plow250" "0.917"
## [1] "plow250" "plow500" "0.847"
## [1] "plow250" "tpi250" "0.755"
## [1] "plow500" "tpi250" "0.708"
## [1] "plow500" "tpi500" "0.783"
## [1] "tpi050" "tpi125" "0.839"
## [1] "tpi050" "tpi250" "0.737"
## [1] "tpi125" "tpi250" "0.93"
## [1] "tpi125" "tpi500" "0.762"
## [1] "tpi250" "tpi500" "0.897"
#
# # distance from Atlantic and False Bay, borderline
# # Ran some models swapping these two variables in and out, and found false bay on average increased R2 more, and overall model R2 are significantly correlated so when one works, either works overall. Based on this, all models use Distance to False Bay.
# cor(topo10[,c('d2at','d2fb')],use='pair') # R = -0.65
#
# # equinox and winter solstice radiation
# cor(topo10[,c('rad080','rad172')],use='pair') # R = 0.95
#
# # all tpi variables with each other across scales, and plow with each other across scales; elevation weakly correlated with tpi and plow
# cor(topo10[,c('elevation','tpi050','tpi125','tpi250','tpi500','plow050','plow125','plow250','plow500')],use='pair') # R >= 0.89
# cor(topo30[,c('elevation','tpi050','tpi125','tpi250','tpi500','plow050','plow125','plow250','plow500')],use='pair') # R >= 0.89
#
How does the Topographic Heat Load index, calculated with different aspect orientations, compare to daily radiation load from a solar model? Maximum values are near the equinoxes, with a minimum at summer solstice and a dip at winter solstice. thl315 has lower correlation as it is furthest from north-facing, so may work for temperature due to hotter afternoon air temp, but is not as strongly correlated with solar radiation per se.
# # How does daily solar radiation compare to topographic heat load
# # correlation max at 0.7 on day 48 (Feb 17) and around day 295 (Oct 21)
# # minimum mid-summer solstice
# head(dlySol)
# head(topo10)
solVthl <- c()
i=1
for (i in 1:366) solVthl[i] <- cor(topo10$thl0,dlySol[,i],use='pair')
plot(1:366,solVthl,type='l')
for (i in 1:366) solVthl[i] <- cor(topo10$thl315,dlySol[,i],use='pair')
lines(1:366,solVthl,col='red')
for (i in 1:366) solVthl[i] <- cor(topo10$thl337,dlySol[,i],use='pair')
lines(1:366,solVthl,col='blue')
The following shows scatterplots of pairwise differences matrices for spatial distance vs. each environmental factor, and the accompanying Mantel test. The plots that are commented out are not significant for spatial correlation, and plots that are significant are illustrated. The random string in the box is an identifier to link this code to the results in the manuscript. Site 32 (Kirstenbosch Race Track) is removed from analysis.
As expected, distance from ocean variables are strongly spatially correlated. This raises questions about their interpretation in the models, but they are generally their for the purpose of capturing spatial gradients, so that could be okay. Elevation and slope are also spatially structured, though the scatterplots suggests less concern. Exposure and hillslope positions are not correlated, indicating we did a good job stratifying across the region.
#
# #--------------------------------------------#
# # SPATIAL VARIATION IN TOPOGRAPHIC VARIABLES #
# # a6STycTJVjqM8vRd8V8h #
# #--------------------------------------------#
library(ade4)
# check matrices are aligned
all(meta$siteID==topo10$siteID)
## [1] TRUE
#head(meta)
sdist <- dist(meta[-32,c('UTM.east','UTM.north')])
hist(sdist,main='Histogram of distances between sites (m)')
selev <- dist(meta$elevation[-32])
plot(sdist,selev)
mantel.rtest(sdist,selev,nrepet=999)
## Warning in is.euclid(m2): Zero distance(s)
## Monte-Carlo test
## Observation: 0.2976591
## Call: mantel.rtest(m1 = sdist, m2 = selev, nrepet = 999)
## Based on 999 replicates
## Simulated p-value: 0.001
sd2at <- dist(topo10$d2at[-32])
plot(sdist,sd2at)
mantel.rtest(sdist,sd2at,nrepet=999)
## Monte-Carlo test
## Observation: 0.4048515
## Call: mantel.rtest(m1 = sdist, m2 = sd2at, nrepet = 999)
## Based on 999 replicates
## Simulated p-value: 0.001
sd2fb <- dist(topo10$d2fb[-32])
plot(sdist,sd2fb)
mantel.rtest(sdist,sd2fb,nrepet=999)
## Monte-Carlo test
## Observation: 0.9611594
## Call: mantel.rtest(m1 = sdist, m2 = sd2fb, nrepet = 999)
## Based on 999 replicates
## Simulated p-value: 0.001
plot(selev,sd2fb)
mantel.rtest(selev,sd2fb,nrepet=999)
## Warning in is.euclid(m1): Zero distance(s)
## Warning in is.euclid(distmat): Zero distance(s)
## Monte-Carlo test
## Observation: 0.2783594
## Call: mantel.rtest(m1 = selev, m2 = sd2fb, nrepet = 999)
## Based on 999 replicates
## Simulated p-value: 0.001
sd2cs <- dist(topo10$d2cs[-32])
plot(sdist,sd2cs)
mantel.rtest(sdist,sd2cs,nrepet=999)
## Monte-Carlo test
## Observation: 0.1044175
## Call: mantel.rtest(m1 = sdist, m2 = sd2cs, nrepet = 999)
## Based on 999 replicates
## Simulated p-value: 0.003
dslope <- dist(topo10$slope[-32])
plot(sdist,dslope)
mantel.rtest(sdist,dslope,nrepet=999)
## Monte-Carlo test
## Observation: 0.0561465
## Call: mantel.rtest(m1 = sdist, m2 = dslope, nrepet = 999)
## Based on 999 replicates
## Simulated p-value: 0.018
dtpi50 <- dist(topo10$tpi050[-32])
#plot(sdist,dtpi50)
mantel.rtest(sdist,dtpi50,nrepet=999)
## Monte-Carlo test
## Observation: 0.02425777
## Call: mantel.rtest(m1 = sdist, m2 = dtpi50, nrepet = 999)
## Based on 999 replicates
## Simulated p-value: 0.176
drad080 <- dist(topo10$rad080[-32])
#plot(sdist,drad080)
mantel.rtest(sdist,drad080,nrepet=999)
## Monte-Carlo test
## Observation: 0.01535659
## Call: mantel.rtest(m1 = sdist, m2 = drad080, nrepet = 999)
## Based on 999 replicates
## Simulated p-value: 0.196
dthl <- dist(topo10$thl0[-32])
#plot(sdist,dthl)
mantel.rtest(sdist,dthl,nrepet=999)
## Monte-Carlo test
## Observation: 0.01342669
## Call: mantel.rtest(m1 = sdist, m2 = dthl, nrepet = 999)
## Based on 999 replicates
## Simulated p-value: 0.26
# #-------------------------------------------------#
# # MAKE 3 CHAR DAYS VAR, for dlySol variable names #
# #-------------------------------------------------#
days <- as.character(1:366)
days[nchar(days)==1] <- paste('00',days[nchar(days)==1],sep='')
days[nchar(days)==2] <- paste('0',days[nchar(days)==2],sep='')
Okay, the modeling strategy is a bit brute force now! First set up four variables that have alternative terms for equivalent topographic variables: radterm = solar radiation; hillterm = topographic hillslope position; regterm = regional position relative to ocean; slterm = slope (only one variable). Plus there’s elevation in every model. Then set up modterms and modnums which contain the names and numbers of each item for all 120 combinations (5 radterm * 8 hillterm * 3 regterm). And finally set up a bunch of lists which will have four items each, holding results for Tmin, Tmax, RHsat.hrs, and VPmax; each item in the list will be a 366 row matrix with model results for each day
#-------------------------------#
# DAILY WEATHER MODELS #
# OHkNt68RXvbTk4jHfczb #
#-------------------------------#
# Run through daily models testing all combinations of elevation, radiation, hillslope, and regional position terms
radterm <- c('dsol','rad080','rad172','rad355','thl315','thl337','thl0')
hillterm <- c('tpi050','tpi125','tpi250','tpi500','plow050','plow125','plow250','plow500')
regterm <- c('d2fb','d2at','d2cs')
slterm <- 'slope'
nModels <- length(radterm)*length(hillterm)*length(regterm)
modterms <- matrix(NA,nModels,5)
modnums <- matrix(NA,nModels,5)
n <- 0
for (i in 1:length(radterm)) for (j in 1:length(hillterm)) for (k in 1:length(regterm)) {
n <- n+1
modterms[n,] <- c('elevation',radterm[i],hillterm[j],regterm[k],slterm)
modnums[n,] <- c(1,i,j,k,1)
}
## running all possible models on daily basis for four dependent parameters: Tmin, Tmax, RHsat.hrs, VPmax
RSQelev <- list() # r-squared, elevation only model
MSEelev <- list() # mean square error, elevation only model
RSQm5x <- list() # highest r-squared, 5 term model
MSEm5x <- list() # mean square error for max r-squared 5 term model
BFm5x <- list() # best-fit model: number from 1:120 indicating which had highest r-squared
RSQm5 <- list() # matrix of 120 r-squared for all 5 term models
MSEm5 <- list() # matrix of 120 mean square errors for all 5 term models
vars <- c('Tmin','Tmax','RHsat.hrs','VPmax')
for (v in 1:length(vars)) {
RSQm5[[v]] <- matrix(NA,366,ncol=nModels)
MSEm5[[v]] <- matrix(NA,366,ncol=nModels)
RSQelev[[v]] <- rep(NA,366)
MSEelev[[v]] <- rep(NA,366)
RSQm5x[[v]] <- rep(NA,366)
MSEm5x[[v]] <- rep(NA,366)
BFm5x[[v]] <- rep(NA,366)
}
Now, run through the four weather variables (v loop) and the 366 days (d loop). Each time through, subset the weather data, and insert the appropriate daily solar radiation from the dlySol data.frame into topo10$dsol. Fit one model with elevation only, and record results. Then fit all 120 possible 5 terms models (i,j,k, loop), saving r-squared and MSE. After running all of them, extract best fit and corresponding r-sq and MSE.
v <- 1
for (v in 1:length(vars)) {
d <- 1
for (d in 1:366) {
#print(c(v,d))
dd <- subset(dw,dw$doy==d)
dd <- dd[match(dd$siteID,meta$siteID),]
dd$yvar <- dd[,vars[v]]
dd$yvar[meta$use4Analyses==0] <- NA
topo10$dsol <- dlySol[,paste('rad_tot_',days[d],sep='')]
# Base model with elevation only
fit1 <- lm(dd$yvar~elevation,data=topo10)
RSQelev[[v]][d] <- summary(fit1)$r.sq
MSEelev[[v]][d] <- sd(fit1$residuals,na.rm=T)
# run all 5 parameter models
n <- 0
i=1;j=1;k=1
for (i in 1:length(radterm)) for (j in 1:length(hillterm)) for (k in 1:length(regterm)) {
n <- n+1
#print(c(d,i,j,k,n))
fit2 <- lm(dd$yvar~topo10[,'elevation']+topo10[,radterm[i]]+topo10[,hillterm[j]]+topo10[,regterm[k]]+topo10[,slterm])
RSQm5[[v]][d,n] <- summary(fit2)$r.sq
MSEm5[[v]][d,n] <- sd(fit2$residuals,na.rm=T)
}
}
## now extract rsq and mse for the best model, and record id # of best model
BFm5x[[v]] <- apply(RSQm5[[v]],1,which.max)
for (d in 1:366) {
RSQm5x[[v]][d] <- RSQm5[[v]][d,BFm5x[[v]][d]]
MSEm5x[[v]][d] <- MSEm5[[v]][d,BFm5x[[v]][d]]
}
}
Now, let’s look at the results! First, how do Rsq for the elevation and the best of the 5-term models vary through the year, for each variable? Interestingly, for all four variables results fluctuate a lot day to day, especially based on elevation alone. Using 5 variables, r-squareds are mostly between 40 and 90%, with lower values for: Tmin in spring and fall, Tmax in summer, RHsat.hrs in winter, and VPmax in summer.
op=par(mfrow=c(2,2))
i=1
for (i in 1:4) {
plot(1:366,RSQelev[[i]],type='l',col='red',ylim=c(0,1),main=vars[i],xlab='DOY',ylab='r-squared')
points(1:366,RSQm5x[[i]],type='l',lwd=2)
}
par(op)
Interestingly, there are only weak correlations between model fit across variables on a given day. The strongest correlation is for the 5-term models for Tmin and Tmax (var 1 vs var 2 in second plot below):
pairs(cbind(RSQelev[[1]],RSQelev[[2]],RSQelev[[3]],RSQelev[[4]]))
pairs(cbind(RSQm5x[[1]],RSQm5x[[2]],RSQm5x[[3]],RSQm5x[[4]]))
My main goal has been trying to make sense of variation in the daily r-squared variables, and next step would be to sort out the meaning of which terms contribute to the best model. Presumably many of the r-squared values are effectively indistinguishable, but that would require quite a lot of examination.
Here are plots of the max r-squared for each variable versus the daily mean for that variable.
Tmin - nothing doing!
plot(dlySummary$Tmin,RSQm5x[[1]],ylab='Max r-squared, Tmin')
Tmax - r-squared are higher on cooler days, i.e. topography is more predictive of spatial variation in Tmax
plot(dlySummary$Tmax,RSQm5x[[2]],ylab='Max r-squared, Tmax')
Tmax variation is also much better explained on days with a low Tmin-Tmax difference, which would presumably be humid and/or cloudy:
plot(dlySummary$Tmax-dlySummary$Tmin,RSQm5x[[2]],ylab='Max r-squared, Tmax',xlab='Tmax-Tmin')
That’s interesting, as I think of the coolest days as often being clear with a potentially high Tmin-Tmax difference. Just a quick check on how those two are related:
plot(dlySummary$Tmax,dlySummary$Tmax-dlySummary$Tmin,xlab='Tmax',ylab='Tmax-Tmin')
plot(dlySummary$Tmin,dlySummary$Tmax,xlab='Tmin',ylab='Tmax')
abline(0,1)
RHsat.hrs - r-squared are higher on days with intermediate number of hours of fog - perhaps these days have the most variation among stations to explain, since clear (=0) and all foggy (=24), there’s not much variation to work with
plot(dlySummary$RHsat.hrs,RSQm5x[[3]],ylab='Max r-squared, RHsat.hrs')
VPmax - not much going on
plot(dlySummary$VPmax,RSQm5x[[4]],ylab='Max r-squared, VPmax')
Inspection of Kirstenbosch weather variables didn’t add much here.
There’s also not much obvious variation among SOM classes from Jasper’s regional weather pattern analysis:
plot(dlySummary$sommin,RSQm5x[[1]],ylab='Max r-squared, Tmin',xlab='Tmin regional SOM')
plot(dlySummary$sommax,RSQm5x[[2]],ylab='Max r-squared, Tmax',xlab='Tmax regional SOM')
# t1 <- table(BFm5x[[1]])
# modterms[as.numeric(names(t1)[which.max(t1)]),]
#
# t1 <- table(BFm5x[[2]])
# modterms[as.numeric(names(t1)[which.max(t1)]),]
#
# t1 <- table(BFm5x[[3]])
# modterms[as.numeric(names(t1)[which.max(t1)]),]
#
# t1 <- table(BFm5x[[4]])
# modterms[as.numeric(names(t1)[which.max(t1)]),]
# modterms[46,]
# modnums[46,]
# plot(modnums[BFm5x[[2]],3])
#
# plot(RSQelev[[1]],RSQm5x[[1]])
# plot(RSQelev[[2]],RSQm5x[[2]])
# plot(RSQelev[[3]],RSQm5x[[3]])
# plot(RSQelev[[4]],RSQm5x[[4]])
#
# plot(MSEelev[[1]],MSEm5x[[1]]);abline(0,1)
# plot(MSEelev[[2]],MSEm5x[[2]]);abline(0,1)
# plot(MSEelev[[3]],MSEm5x[[3]]);abline(0,1)
# plot(MSEelev[[4]],MSEm5x[[4]]);abline(0,1)
#
# head(dlySummary)
# boxplot(MSEm5x[[2]]~dlySummary$sommax)
# boxplot(MSEm5x[[1]]~dlySummary$sommin)
#
# plot(RSQm5x[[1]]~dlySummary$Kps.mean_hPa)
# plot(RSQm5x[[2]]~dlySummary$Kps.mean_hPa)
# plot(RSQm5x[[3]]~dlySummary$Kps.mean_hPa)
# plot(RSQm5x[[4]]~dlySummary$Kps.mean_hPa)
#
# plot(RSQm5x[[1]]~dlySummary$Kwspd.mean_m.s)
# plot(RSQm5x[[2]]~dlySummary$Kwspd.mean_m.s)
# plot(RSQm5x[[3]]~dlySummary$Kwspd.mean_m.s)
# plot(RSQm5x[[4]]~dlySummary$Kwspd.mean_m.s)
#
# plot(RSQm5x[[1]]~dlySummary$Tmin)
# plot(RSQm5x[[2]]~dlySummary$Tmax)
# plot(RSQm5x[[3]]~dlySummary$RHsat.hrs)
# plot(RSQm5x[[4]]~dlySummary$VPmax)
#
# plot(RSQm5x[[1]]~I(dlySummary$Tmax - dlySummary$Tmin))
# plot(RSQm5x[[2]]~I(dlySummary$Tmax - dlySummary$Tmin))
# plot(RSQm5x[[3]]~I(dlySummary$Tmax - dlySummary$Tmin))
# plot(RSQm5x[[4]]~I(dlySummary$Tmax - dlySummary$Tmin))
#
# plot(I(dlySummary$Tmax - dlySummary$Tmin)~dlySummary$Tmin)
# plot(I(dlySummary$Tmax - dlySummary$Tmin)~dlySummary$Tmax)
# plot(I(dlySummary$Tmax - dlySummary$Tmin)~dlySummary$RHsat.hrs)
# plot(I(dlySummary$Tmax - dlySummary$Tmin)~dlySummary$VPmax)